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Abstract 

We address a linear fractional differential equation and develop effective solution methods 
using algorithms for inversion of triangular Toeplitz matrices and the recently proposed 
QTT format. The inverses of such matrices can be computed by the divide and conquer 
and modified Bini's algorithms, for which we present the versions with the QTT approx- 
imation. We also present an efficient formula for the shift of vectors given in QTT format, 
which is used in the divide and conquer algorithm. As the result, we reduce the com- 
plexity of inversion from the fast Fourier level (n log n) to the speed of superfast Fourier 
transform, i.e., O(log 2 n). The results of the paper are illustrated by numerical examples. 
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1. Introduction 

Equations involving derivatives of fractional order are of great importance, due to their 
role in mathematical models applied in mechanics, biochemistry, electrical engineering, 
medicine, etc., see [|TTll8llT8l. In this paper we present a superfast algorithm for the nu- 
merial solution of the linear equation 

D^(t)=F(t,y(t))=my(t)+f(t), < t < T, y(0)=y o , (1) 

where < cx < 1 is the order of the fractional operator, m e K is a constant referred to 
as mass, and f (t) is a sufficiently well-behaved forcing term. For cx = 1/2 this equation is 
a scalar version of the Bagley-Torvik equation O, which is used in the modelling of vis- 
coelastic materials. The definitions of Caputo derivative D * can be found in many sources, 
e.g. flT2~ll40|] , and are presented in appendix for the convenience. 

The classical result of Diethelm | |12t Lem. 6.2] allows us to rewrite ([T]) in the form 

1 ft (t-s) a - 1 (mu(s) + f(s))ds, (2) 



y(tj = y + 



o 
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where F(<x) = e _t t a_1 dt is the gamma function. Eq. Q is the weakly singular convolu- 
tional Volterra equation of the second kind with the Abel-type kernel. Volterra equations 
of second kind are well-studied and are proven to have a unique continuous solution for 
< t < T, see, e.g. ||26l Thm. 3.2]. The solution is asymptotically stable if m < (see ||21||) 
which we will always assume in this paper. 

For certain forcing terms, the solution of §2§ can be found using series methods. In 
a general framework, we can discretize ((2J using a collocation or Galerkin method and 
numerically solve the resulted linear system. This matrix approach to fractional calculus 
was brilliantly presented by I. Podlubny in ||41|. In this paper we consider the collocation 
method and assume that y (t) is approximated by a piecewise-linear function on a uniform 
grid tj = jh, j = 0, . . . , n, where h = T/n. The stability of collocation methods for fractional 
equations was studied in B6J.0 and an error analysis can be found in |[T4|. The discretized 
equation is the following 

yj =yo + -p^^w j)k (ray k + f k ), j = l,...,n, 
1 [0L> k=0 

where yj = y (tj), f k = f(t k ) and Wj )k are quadrature weights, defined by integration of 
piecewise-linear basis functions with Abel-type kernel, i.e., 
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(j-1) a+1 -(j-a-1)j a , k = 0, 

Wj, k = - { 0-k-1) a+1 -2(j-k) K+1 + 0-k+1) a+1 , 1<k<j, 
a(a + 1) ] ^ k = j> 

Finally, we obtain the linear system Ay = b with triangular Toeplitz matrix and the right- 
hand side defined as follows, 



kyk = b h j = 1,...,n, (3) 



k=1 



n = j 1 - T m > P = °) 

v \ -ym((p-1) a+1 -2p a+1 + (p + 1) a+1 ), p>0, 

bj =y + y \ Y_ w j|k f k + Wj i0 (myo + f ) J , 

where y = h a /r(oc + 2). 

The numerical scheme we use is analogous to the fractional Adams method proposed 
in (HI for a general (e.g. nonlinear) function F(t,y(t)). The method is developed as a gen- 
eralization of the Adams-Bashforth-Moulton scheme from the classical numerical analy- 
sis of ordinary differential equations and a detailed error analysis is provided. The com- 
plexity of the fractional Adams method in the nonlinear case is 0(n 2 ). To reduce this com- 
plexity, we can take into account the decay speed of the Abel kernel k(s) = s a 1 of the 
integral in The so-called fixed memory principle [40ll39l and more accurate nested mesh 
method | |T7l IT3l are based on truncation and approximation of the tail of the integral 
respectively, and have almost linear complexity w.r.t. n. We revise these methods in Sec.|2j 

For linear F(t,y (t)), the problem writes as the linear system ([3]), which can be solved 
using well-developed algorithms for the inversion of triangular Toeplitz matrices, or tri- 
angular strip matrices, as they are referred in [ j4"T1 ]. These methods are recalled in Sec. [3j 
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and have the asymptotic complexity of the fast Fourier transform (FFT) algorithm, which 
is O(nlogn). 

Recently, a superfast Fourier transform algorithm was proposed in ||15|, based on the 
approximation of vectors in the quantized tensor train (QTT) format P5ll2lf. The method 
can be considered as a classical model of quantum superfast Fourier transform algo- 
rithm [ 116] . and has a square-logarithmical complexity O(log 2 n) for a certain class of vec- 
tors, for which such a model is efficient. This class of vectors is partially established in | |42 || 
and include, for example, vectors with sparse Fourier image. The numerical experiments 
provided in Sec. [5] show that the Abel kernel t(s) = s 1_a is efficiently approximated by 
the QTT format for all < <x < 1 with accuracy up to the machine threshold. Based on 
this observation, we propose the superfast inversion algorithm for the triangular Toeplitz 
matrix ([3]), using the QTT approximation. 

The numerical experiments provided in Sec. |5]justify the accuracy and sublinear com- 
plexity of the method proposed. 



2. Numerical method with logarithmic memory 



In [ j40] and [39| the author describes an approach to the numerical integration involved 
in solving a fractional problem whereby the first part (or tail) of the integral is ignored (i.e. 
assuming the value of the integral over this region is negligible) and so the memory of the 
system is truncated at some point. The error introduced via this process is described in 
[40J for Riemann-Liouville fractional derivatives. In | |17 || the authors consider the error 
that is introduced when this approach is applied to problems expressed with respect to 
the Caputo fractional derivative. The authors show that by introducing a finite memory 
of fixed length T for the Caputo derivative we introduce an error of the form 
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rt-T 



rn 



cx 



y'ls) 



rdS 



(4) 



Letting sup se[0>t] \y'{s)\ = M then B 



E < 



M (t 1 - tx -T 1 - K ) 
r (2 - a) 



(5) 



So for a fixed memory T < t we have a loss of order such that the error does not tend 
towards zero as the stepsize approaches zero. Indeed, the authors in | |17 || highlight that in 
order to preserve the order of the method we would need to choose T so that (for a fixed 
error bound F) we have 

which introduces a computational cost — precisely what the fixed memory principle is 
trying to avoid. To overcome this it is proposed in ||17||, and described further in [13J, 
that the fixed memory principle is amended so that the region of integration [0, t] is de- 
composed into a sequence of finite-length intervals with differing stepsizes. So as we 
move 'backwards' along the interval from t to the subintervals use coarser and coarser 
stepsizes, except possibly for some small sub-interval near zero due to the length of this 
subinterval not being an exact multiple of the current stepsize — in such circumstances 
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Figure 1: Example of the grids used at subsequent steps of the nested mesh method (from top to bottom). 
On each line the active points of the grid are shown by black and non-active by grey dots. 

the authors suggest a couple of alternative approaches for this subinterval, one such al- 
ternative being the use of the original stepsize. Such a nested mesh approach (see actual 
mesh on Fig. [T]) is possible due to the scaling properties of the fractional integral, which 
are discussed in | |13 || and ||17|. Thus the weights (of the Adams-type method described 
earlier) for calculating Q£f (nh) « I a f (nh) with a stepsize H can be used to calculate 
n^ Ph f (cu p nh) « I a f (ntu p h) using a stepsize of cu p h. The authors [17] define, for h G IR + , 
the mesh M h by M h = {hn, n e N}. If co, r, p e N, w > 0, r > p, then M w r h c M.^. The 
authors then decompose the interval [0, t], for fixed T > in the following way: 

[0, t] = [0, t - w m T] U [t - cu m T, t - w m - ] T] U • • • U [t — cuT, t — T] U [t — T, t] (7) 

where m e N is the smallest integer such that t < cu m+1 T. A step length of H is used 
over the most recent time interval [t — T, t] with successively larger step sizes over earlier 
intervals, as follows. Let t, T, h e E, uj m+1 T > t > cu m T, t > 1, h > with t = nh for some 
n E N. The integral can be rewritten as 

m-1 

if 0)t] f(t) = i^f (t) + Y_ ift-o^tWW + nf W (8) 

i=0 
m-1 

= l^ t] f(t) + Y_ ^IgUou^nf i^t) + ^ll t _ wmJ] f (cu-t) , (9) 

i=0 

where I^_ at _ b] f(t) = ^ ^yr=^ ds. The authors also show the following: 

Theorem 2.1. /II 71/ The nested mesh scheme preserves the order of the underlying rule on which it 
is based. 

In addition, whilst the computational cost of the full-memory approach is of order 0(N 2 ), 
the nested-mesh approach has order 0(1M log N). 



3. Inversion of triangular Toeplitz matrices 



3.1. Basic properties of triangular Toeplitz matrices 

Let T n be a set of lower triangular Toeplitz n x n matrice^J i.e., 

A e T n & A=[a(j,k)]^ , a(j,k) = a(j -k), a(p) = 0, p < 0. 

It is easy to check the following properties of T n . 

LAG T n , B G T n AB G T n ; 

2. A G T„, B G T n 4 AB = BA; 

3. A G T n , a ^ 4 A" 1 G T n . 

By the last property, the inverse matrix B = A -1 , as well as all matrices from T n , is de- 
fined by its first column. The standard solution method for triangular linear systems has 
complexity 0(u 2 ) and yields the following trivial formula 



M0) = -jL b(j)=-J-^b(j-k)a(k), j = l,...,n-l. (10) 
a(0) a(0) ^ 



For A, B G T n , the product X = AB G 7 n and is also defined by the first column x = Ab. 
Therefore, matrix-by-matrix multiplication in T n is equivalent to the multiplication of a 
vector by the Toeplitz matrix, i.e., discrete convolution x{)) = Y-k=o a 9 ~~ k)b(k). A naive 
computation by this formula requires (u 2 ) operations, but it is well-known that it can be 
computed in Ofulogn) operations using the fast Fourier transform (FFT) algorithm | |19t 
I10 ||. To recall this, we note that each n x n Toeplitz matrix T is the leading submatrix of 
some 2n x 2n circulant matrix 



C 



T * 
* T 



C = [c(j, k)] , where c(j, k) = c(j — k mod 2n), 



and all circulant matrices are diagonalized by unitary Fourier matrix as follows (see 
cf. B2Q1) 

C = F*AF, A = Vlndiag{¥c). 

Therefore, multiplication by C and hence by T can be performed by 3 FFTs of size n with 
complexity (u log n) . 

The inversion of triangular Toeplitz matrices has asymptotically the same complexity, 
i.e., cM(n), where M(n) denotes the complexity of matrix multiplication. The modern 
highly-improved inversion algorithms reduce the constant to the level from c = 1 .4 to 
c = 1 .5, see, e.g. ||30||. We now recall the classical algorithms, which have slightly larger 
constant c, but are much more simple and easy to follow. In Sec.|4]we will adjust the classi- 
cal inversion algorithms to use the compressed format for the approximate representation 
of matrix, reducing the complexity to sublinear w.r.t. n. 



2 Here and further we write matrix and vector indices in round brackets instead of putting them as sub- 
scripts, in order to introduce the convenient notation for QTT representation later. 
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3.2. Divide and conquer method 

To benefit from the Toeplitz structure and reach O(nlogu) complexity for the inver- 
sion algorithm, we can use the divide-and-conquer strategy. This was noted in [29] and 
developed in J9J. It is easy to check that if 2n x 2n lower triangular Toeplitz matrix A' G T 2n 
is partitioned to n x n matrices, the inverse matrix writes as follows 



A' 



A 

C A 



fA'l 





" A- 1 






"A- 1 






I 






-A" 1 CA- 1 


A- 1 






A- 1 




-CA" 1 


I 



(11) 



2 d and A d G T 2 d, this formula 



where A G T n , A 1 G T n and C is a Toeplitz matrix. If n 

yields the reccurent method to compute A^ 1 . We start from some small do and use ( fToj ) to 



compute the inverse of 2 d ° x 2 d ° leading submatrix Aa G T 2 d . Then we subsequently ap- 
ply ( [TT) and compute A d 1 in (d— d ) steps. Each step requires to compute the first column 
of AfXtA^ 1 with 2 t x 2 t Toeplitz matrices A^ 1 and C t , where t = d + 1 , . . . , d. Each mul- 
tiplication is done in 0(t2 t ) operations, which summarizes to 0(d2 d ) = O(nlogn) overall 
complexity. More accurately, the cost of the divide and conquer algorithm is smaller than 
1 2 FFTs of size n. 



3.3. Bini's and related approximate methods 

In order to reduce the number of FFTs used in computations and obtain algorithm 
with better parallel performance, the approximate method to compute A 1 for A G T n was 
proposed in [5J. It is noted that T n is the algebra generated by the matrix H G T n with unit 
elements on the subdiagonal and zeros elsewhere, i.e., transposed Jordan block with zero 
diagonal. Therefore, A G T n with first column a = [a(j)]^T 1 is written as A = Y.^=o i(j)H'. 
The idea is to add a small element £ n at the top right corner of the matrix and substitute 
H by H £ = H + e n eje n _ 1 . It is easy to check that Dc^D^ 1 = eC, where D £ = diag{e 5 '}jr o \ 
and C = Hi = H + ejen i generates the algebra of circulant n x n matrices. Then A and 
A -1 are approximated as follows 

n-1 /n-1 \ 

A « A £ = Y_ a(j)H* = D7 1 X a(j)£ j q D £ = D^QD, = D f T 1 F*A £ FD e , 

j=o \j=o / ( 12 ) 

A- 1 « A- 1 = D7 1 F*A7 1 FD £ , where A £ = v^diag(Fa £ ), a e (j) = a(])e\ 

The first column of A7 1 is computed using two FFTs of size n. 

This idea was revised in |25|, where it was proposed to apply Bini's algorithm to the 
first column a of matrix A padded with n zero elements. The revised version of Bini's 
algorithm requires two FFTs of size 2n and has better accuracy properties. 

3.4. Newton iteration 

The classical Newton iteration 



B k+1 = 2B k - B k AB k (13) 

was proposed in | |45 || for the computing the inverse A 1 of a nonsingular matrix A. It con- 
verges quadratically if initial guess B is s.t. ||I — AB || < 1 in any operator norm of a 
matrix. In JU it is shown that for B = (jA* with some small real \i Newton iteration con- 
verges to the inverse A 1 of a nonsingular or pseudoinverse A^ of a singlular matrix A. 
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In [ j37fl further deep analysis is provided, for instance, it is shown that pr 1 = HA^HAH^, 
is a good and reliable choice. In relation to Toeplitz and related structured matrices, the 
Newton iteration with approximation of the result on each step was developed using the 
concept of displacement ranks | |38 || and tensor product approximations |[2"2~ll3"Tll36l] . 

For A G T n the choice of initial guess B = fxA* is not effective, since A* ^ T n and we 
can not perform iterations with G T n , which grants low storage and fast multiplication. 

If B G T n , every Newton iteration costs two convolutions, i.e., 6 FFTs of size 2n. 

Remark 1. A single Newton iteration for lower triangular Toeplitz matrices is slower than the 
divide and conquer method. 

It is not easy to provide a good initial guess B G 7 n for which the Newton iteration with 
a given matrix A G T n converges in one or few steps. However, Newton iteration can be 
used to improve the accuracy of matrix B « A 1 , B G 7 n computed by other means, if 
|| I — AB || < 1 . For instance, we can note the following relation between the divide and 
conquer method and Newton iteration. 

Remark 2. For matrix A' G 7i n defined in ( |TT| , the Newton iteration ( [13] ) with initial guess 



Br 



A- 1 0' 




A G T n , 



gives B i = ( A') 1 , i.e., converges in one step and is equivalent to the divide and conquer method ( |TT] >. 

Therefore, for A G T n , divide and conquer method is always better than the Newton iter- 
ation, which reduces to the divide and conquer method in the special case. 

3.5. Decay of the elements of inverse matrix 

It is instructive to look at the decay profiles of the elements of a triangular Toeplitz 
matrix Q and its inverse, see Fig. [2] There is a jump in magnitude between diagonal and 
subdiagonal elements, i.e., 

a 1 - (ym)-' h a 



cn 2 a+1 -2 ' 7 r(a + 2)' 

where the numerator increases when n — > oo, h, — > and tends to one when m — > — oo. 
After the jump, elements decay polynomially, i.e., a p ~ p a 1 for p > 1 . For the inverse ma- 
trix the behaviour is the same for a certain (possibly very long) set of elements. However, 
after certain point the rate of decay changes from 1 — a to 1 + a, i.e., b p ~ p a 1 for p > P. 
The bend point P which is obtained from the experiment, is the monotonically decreasing 
function P = P(yrri), i.e. the larger is the initial jump, the later the decay of element of 
the inverse matrix switches to faster rate. The observed behaviour of elements of inverse 
matrix allows us to predict the upper bound for the norm of the second half of vector, us- 
ing the information about the first half. We will use this property in the next subsection, 
where the divide and conquer algorithm will be adapted for the vectors approximately 
given in the low-parametrical tensor-structured format. 
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Figure 2: Decay profiles of triangular Toeplitz matrix (|3) (left) and its inverse (right) for n = 2 28 , h = T/n, 
T = 10 and <x = 0.5 (top) and a = 0.8 (bottom) and for different mass m. 



4. Inversion of triangular Toeplitz matrices using QTT approximation 

4.1. Tensor train and quantized tensor train formats 
A tensor is an array with d indices (or modes) 



[a(ki,...,Tcd)], 



The tensor train (TT) format 



kp — 0, . . . , Tip 1 , 
for the tensor A reada^ 



l,...,d. 



a(k l ,k 2 ,...,k d )=A[ 1 X 2 2 ) ...A (d: 



k d > 



where each a{, p ' is an r p _i x r p matrix. Usually the border conditions To = r d = 1 are imposed 
to make every entry a(k| , . . . , k d ) a scalar. However, larger r and r d can be considered and 
every entry of a tensor A = [a(k| , . . . , k d )] becomes an r x r d matrix. Values r , . . . , r d _i 
are referred to as TT-ranks and characterize the separation properties of the tensor A. Three- 
dimensional arrays A' p ' = [A{£ ] are referred to as TT-cores. 



(14) 



3 We will often write the equations in elementwise form, which assumes that all indices run through all 
possible values. 
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Figure 3: Effective QTT rank of vector [k"" 1 ] (top left), vector [(k - 1 ) a+1 - 2k a+1 + (k + 1 ) a+1 ] (top right), 
first column of matrix <j3j (bottom left) and its inverse (bottom right) w.r.t. parameter a and relative approx- 
imation accuracy e in Frobenius norm. Problem size n = 2 28 , maximum time T = 10, mass m = — 10 +s 



To apply the TT compression to low dimensional data, the idea of quantization was 



proposed fl33l l2"4"l . We will explain the idea for a one-dimensional vector a = [a(k)]£ =( j, 
restricting the discussion to n = 2 d . Define the binary notation of index k as follows 

d 

k = k 1 ...k d = ^kp2 p " 1 > kp = 0,1. (15) 
The isomorphic mapping k <-> (ki , . . . k d ) allows us to reshape a vector a = [a(k)] into the 



d-tensor A = [d(ki , . . . , k d )] . The TT format ( 14 1 for the latter is called the QTT format and 
reads 

a(k) = a(k!...k d ) = d(k! , . . . , k d ) = A™ . . . A^. (16) 

This idea appears in [33| in the context of matrix approximation. In | |24 || the TT format ap- 
plied after the quantization of indices was called the QTT format and applied to a class of 
functions discretized on uniform grids, revealing the impressive approximation proper- 
ties. In particular, it was proven that the QTT-ranks of exp x, sin x, cos x, x p are uniformly 
bounded w.r.t. the grid size. For the functions e ax , x a , etc., similar properties 

were found experimentally. 

The QTT separation function of the function x a 1 discretized on a uniform grid, is 
particularly important for us, because it motivates the use of the QTT approximation to 
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Figure 4: Runtimes of TT-SVD and TT-ACA algorithms for the approximation of matrix (|3j in the QTT 
format w.r.t. size n. (left) a = 0.1, (right) a = 0.9. 



develop the superfast algorithms for the solution of fractional differential equations. In the 
numerical experiment we found out that the QTT ranks are very moderate for all < oc < 1 
and for accuracy up to the machine threshold. The same holds for the first column of 
matrix Q as well as for its inverse. On Fig. [3] we show the effective (average) QTT rank 
w.r.t. e and ex. We note that the effective rank does not overcome 1 0, even for very large 
grids up to n = 2 28 . 

To construct a superfast algorithms in the QTT format we first have to compress the 
data to this format using the algorithm with the sublinear complexity. The original TT- 
SVD algorithm proposed in [35] requires all elements of tensor and therefore does not 
suit for this purpose. To compress matrix Q to QTT format, we apply cross interpolation 
algorithm TT-ACA proposed in [43 j. This method computes the approximation using 
only a few elements of the original array, and does not require all elements to be computed. 
The comparison of runtimes of TT-SVD and TT-ACA algorthms is given on Fig. |4j The 
time that is required to choose the good subset of elements for the interpolation depends 
on the structure of data, which is defined by parameters cx and m. It is easy to see that the 
behaviour of data is less regular for the large cx and mass, which leads to larger runtimes of 
TT-ACA. Nevertheless, we clearly see that TT-ACA outperforms TT-SVD for all examples 
and has sublinear complexity w.r.t. problem size n. 

4.2. Fourier transform and convolution in QTT format 

Inversion algorithms for triangular Toeplitz matrices A 6 T n recalled in Sec.[3]are based 
on two main operations: Fourier transform and discrete convolution. The radix-2 rec- 
curent relation which was known to Gauss | |19| and lays behind the famous Cooley-Tuckey 
FFT algorithm [10] perfectly matches the multilevel structure of QTT format, resulting in 
the QTT-FFT algorithm | jT5| . For a vector of size n = 2 d given approximately in QTT for- 
mat ( 16 ), the QTT-FFT computes the Fourier transform with complexity ( d 2 R 3 ) , where R 
is the maximum QTT rank of the input vector, Fourier image, and all intermediate vectors 
of the algorithm. 

The discrete convolution, i.e., multiplication by Toeplitz matrix, can be performed by 
three Fourier transforms with complexity ( d 2 R 3 ) , where R bounds the QTT ranks of both 
vectors to convolve as well as their Fourier images. As shown in (23]], the convolution 
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c = a * b of two vectors with QTT ranks bounded by r a and Tb , can be written in QTT form 
with QTT ranks bounded by 2r a r b . This representation has large QTT ranks, which can be 
reduced to the value bounded by r c < 2r a r b using some TT-truncation algorithm. We can 
use SVD-based algorithm proposed in | |35 || or iterative DMRG-type approach proposed 
in ||34||, resulting in convolution algorithms with O(dr^r^) and 0(d(r Q + r b + r c )r Q r b r c ) 
complexity, respectively If r Q « r b ~ r c « R, the QTT-FFT and DMRG-based convolution 
algorithms have complexity 0(d 2 R 3 ) and 0(dR 4 ), respectively Therefore, we can not say 
in general which approach is better, even in the simplest case of almost equal QTT ranks. 
This will be established in numerical experiments. 

4.3. Shifts of vectors in QTT format 

The convolution algorithm proposed in [23J is based on the remarkable property of 
shift matrices L 6 T 2 d and U = L T , where the first column of L is I = (0, 1 , 0, 0, ... , 0) T . It 
is shown in [23J that all matrices L p , p = 1 , . . . , 2 d — 1 have all QTT ranks two. Hence, if a 
vector a has QTT ranks ri , . . . , i , then the right shifted vector b = L p a for all p has QTT 
ranks not larger than 2x\ , . . . , 2x&-\ . The same holds for left shifts c = U p a. 

In the following theorem we improve this result for vectors, shifted by one element. 



Theorem 4.1. Let a = [a(k)] k=0 has the QTT representation ( [16] ), then the vector 

b=[x a(0) ... a(2 d -2)] T 



has the QTT representation b(k) = b(kj 



B[ d) with the following cores 



B = 



bS 1] = 



A 



B 



B 



(p) 



(p) 



(p) 



A 



(p) 



B 



B 



(d) 

o — 


A o 


X 


l d) = 


'A[ d] '' 







(17) 



where p = 2, . . . , d — 1 and b q = A\'' ... A^ 13 A^' for q = 2, . . . , d. Similarly, the vector 



= [ail 



a(2 d 



1 



has the QTT representation c(k) = c(ki . . . k d ) = C^ 1 . . . with the following cores 



r' 1 ) _ 



C 



A 



r (p) _ 






r (d) 
u — 


A 


. C P . 




Cd 


p(p) _ 


A?' - 
1 


> 


r (d) 
H — 


A< d) " 








. V . 



(18) 



where? = 2,..., d-1 and c q = A^ IJ . . . A^ A\ q > for q =2,...,d. 



Proof. We check (17 \ straightforwardly. For k = it holds 



b(0)=B 



B 



(d) 



A 



(p) 



A 



= X 
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For k = ki k 2 . . . k d with ki = 1 it holds 



b(k) = b(1k 2 ...k d ) 



B (1) B (2) B (d) 
1 k2 kd 



A 



A 



(2) 



k 2 
* * 



A 



(dV 



kd 

* 



A m A (2) A (d) 



a(0k 2 ...k d ) = a(k-l), 



where "*" denotes arbitrarily zero or non-zero element. For k = kik 2 k 3 . . . k d with ki = 
and k 2 = 1 it holds 



b(k) = b(01k 3 



(3] 



,k„) =B "Br J B 



= b 2 A k3 



,A 



(d) 



A l A A k 3 



B k d - 

(d) 



1 




"aS 2) ■ 




r A 0) - 








M 




* * 




* 



Finally, for k = k] k 2 k 3 . . . k d with ki 

b(k) = b(^_g ik p+1 ...k d ) = 

p— 1 zeros 



A k u ; = a(10k 3 ...k d ) = a(k-1) 
= kp_i = and kp = 1 it holds 



Bi 1] Bi 2) ...Br i) B^Br i) ...B[ d) 



1 




\a (1] ' 
A o 




r A (p-i) - 
A o 




"aS p] i 




A (P + i) 








1 




1 




.b p _ 




* * 




* 



A 



(di 



A 



a(l_^J0k p+1 

p— 1 ones 



•• A l A A k p+1 

a(k-1). 



A 



(d) 



Equation (18 1 is verified in the same way 

4.4. Divide and conquer algorithm in QTT format 

We are now ready to present the version of divide and conquer algorithm which op- 
erates with data given approximately in QTT format. Let n = 2 d and consider A G T n 
defined by the first column a(k) which is represented in the QTT format (16 1. As previ- 
ously let A t denote 2 l x 1} leading submatrix of A. For small do we can invert A do using 
standard divide and conquer method and approximate the first column of A^ 1 in QTT 
format using SVD-based algorithm [35| . 

Now suppose that for some t the first column of A^ 1 is computed in QTT format, and 
we have to compute the QTT approximation of the first column of A t ] ] 1 1 , using the recur- 
sion (11 ). It is necessary to describe the Toeplitz matrix C which lies in the lower part of 
A t+1 . The first column of C is c + = [a(2 t ), a(2 t + 1 ), . . . , a(2 t+1 — 1 )] T and has the following 
QTT representation 



c+fa • • • let) = a(ki . . . k t + 2*) = AjJ J Ag . . . aJJa^A 



(t+2) 



...A 



(d) 



(19) 



The first row of C is c_ = [a(2 t ), a(2 t — 1 ), . . . , a(l )] T . To construct the QTT representation 
for c_, first write the QTT representation for a = [a(0), . . . , a(2 t — 1 )] T , which is 



A 0) A (2) a' 1, A 
A k, A k 2 • • • A k t A 



(t) a (t+1) 



o 



...A 



(d] 



Then apply the 'pull' operation and construct QTT format for a' 
follows 



[a(1 ),..., a(2 t )] T as 



a'(ki-..1ct) 



r (1) r (2) r (t) , (t+1) , (t+2) ,[d) 
^k, ^k 2 • • • C k t A A • • • A > 
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Algorithm 1 Divide and conquer in QTT format 



Require: A G T n , n = 2 d , given by vector [a(k)]£~(j in QTT format ( fT6| ) 
Ensure: B = A -1 G T n given in QTT format 

1: For small do, compute the first column of 2 d ° x 2 d ° leading submatrix A do . Compute 
B do = A^ 1 by ( [T0] > and approximate in in QTT format by TT-SVD algorithm ||35|. 

2: for t = d , . . . , d — 1 do 

3: Compute the first row and column of matrix C t in (111 in QTT format by (19) 
and ((20]). 

4: 

5: Combine the first column of B t and first column of B t C t B t given in QTT format as 
follows 

b(kTT^) = B£ . . . B« g(kT^) = . . . 
to the single vector b' in QTT format, which is defined as follows 



Compute the first column of B t C t B t by two convolutions in QTT format, see Sec. 4.2 



b / (k 1 ...k t k t+1 ; 



B 



J k, 



B 



(2) 
k 2 



] k 2 



B 



r (t] 
b k t 



1 - k t+1 
kt+i 



6: Apply TT-truncate algorithm to V to reduce the ranks of QTT representation. 
7: end for 



where TT-cores Q. are defined by ( [T8] ). Finally, revert the ordering of elements in the 
vector a' to obtain the QTT format for c as follows 

MkTT^ko = c\\ cS\ . . . cj VT'aT 5 • • • < d) - ( 2 0) 

We summarize the above steps in Alg. [T| Note that the workhorse of divide and con- 
quer method is the discrete convolution in QTT format, which can be performed by two 
different methods. This results in two variants of algorithms with different performance, 
which will be studied in numerical experiments. 

4.5. Modified Bini's algorithm in QTT format 

The implementation of ( [12] ) in the QTT format is very straightforward. It is enough to 
mention that the QTT format of vector [e']]^ 1 , n = 2 d , has QTT-ranks one (see B2H), since 

£ j _ gjm-jd _ £ ii £ 2 )2 £ 2d ~'jd 

Therefore, multiplication of a vector in QTT format by diagonal matrix D £ requires only 
the appropriate scaling of TT-cores. By Alg. [2] we present the QTT version of the modified 
Bini's algorithm J25J. Alg. 2]. The algorithm includes two Fourier transforms in the QTT 
format which can not be substituted by discrete convolution. Note that this algorithm 
contains two approximation errors: 

• The first comes form original approximation of triangular Toeplitz matrix A by the 
diagonally scaled circulant matrix A £ . The accuracy of this approximation is gov- 
erned by parameter e. According to the numerical tests made by the authors of [25J, 
the good choice for Bini's and modified Bini's methods are e n = 0.5 x 10 8 and 
e n = 10~ 5 , respectively. 
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Algorithm 2 Modified Bird's method in QTT format 

Require: A 6 T n , n = 2 d , given by vector [a(k)]^ in QTT format ( |T6| ) 
Ensure: B £ = A7 1 ps A -1 G T n given in QTT format 
1: Choose < e < 1 and let ft(lc) = £ k a(k) for k = 0,...,n- 1 and d(k) = for 
k = n, . . . , 2n — 1 . The QTT representation of 0. is the following 

a(k) = a(k 1 ...k d k d+1 ) = Ag 1 . . . A< d d ' (i - k d+1 ), AjJ = aJJ, p = i , . . . , d. 



2: Apply QTT-FFT | |15| to compute the size-2n Fourier transform A = v^nFQ. 

3: Apply Newton iteration ( fT3] ) to compute c = A 1 in the QTT format. Each iteration 

step includes the pointwise (Hadamard) multiplication of vectors in QTT format and 

TT-truncation to reduce the QTT-ranks. 
4: Apply QTT-FFT again to compute the size-2n Fourier transform b = F*c/v / 2n in the 

QTT format 

6(k)=6(k,...k d k a+ ,)=ef , , , ...e« , e£::'- 

5: The QTT representation of the first column of B £ is the following 

b £ (k) = b £ (kT^) = b£> . . . B[ d X d+1) , ^ = r^'^w P = 1 , . . . , d. 



• The second error comes from TT-truncation algorithm applied in the QTT-FFT al- 
gorithm and on each step of Newton iteration. The threshold parameter of TT- 
truncation should be usually smaller than e n in order to maintain the accuracy of 
the result after diagonal scaling. 



5. Numerical experiments 

5.1. Timings of inversion algorithms 

On Fig. [5] we show the runtime of inversion algorithms for triangular Toeplitz matri- 
ces in full and in the QTT format w.r.t. problem size and for different parameters <x, m. 
Standard inversion algorithms have the O(ulogn) complexity which depends only on 
problem size. Quite contrarily, the complexity and runtime of QTT algorithms depend on 
QTT-ranks of input and intermediate vectors, which are sensitive to the fractional order 
a, mass m and step size H. They also depend crucially on the method used to compute the 
discrete convolution in the QTT format. We can note that the divide and conquer algo- 
rithm [l] which uses QTT-conv algorithm [ j23l is always significantly faster than the same 
method which uses QTT-FFT algorithm [15] to compute the convolution. However, QTT- 
FFT works well in modified Bini's algorithm |2j which appears to be the fastest method 
when mass is small in modulus. For large mass the divide and conquer algorithm [T] with 
QTT-conv is preferable to the modified Bini's algorithm |2j For mass m 1 these meth- 
ods have the same asymptotical complexity. 

From Fig. [5]it can be easily seen that the QTT algorithms are asymptotically faster than 
the algorithms in full format. For practical computations it is very important at which size 
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log ]0 (time, sec.) 



log, „ (time, sec.) 




Figure 5: Runtimes of divide and conquer algorithm (solid lines) and modified Bini's algorithm (dashed 
lines) for the inversion of triangular Toeplitz matrix (|3j in full and in the QTT formats (grey and black lines, 
respectively) w.r.t. problem size n and step size h = T/n. Fixed maximum time T = 10, fractional order 
a = 0.2 (left) and a = 0.8 (right), mass m = -10~ 5 (top), m = -10° (middle), m = -10 5 (bottom). 
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log, n 



Figure 6: Accuracy of the solution of the test problem 1 21 1 in the relative Frobenius norm w.r.t. problem size 
n and for different fractional parameters a. Fixed maximum time T = 10 (left) and T = 10 5 (right). Mass 
m = — 1. 



n there is a crossover point, i.e., the minimum value of n for which the QTT algorithms are 
actually faster than the algorithms in full format. Numerical experiments show that for a 
wide range of parameters <x and m the crossover point between full and QTT divide and 
conquer methods is log 2 n ~ 20. This value is about the same as the crossover point be- 
tween FFT and QTT-FFT algorithm applied to signals with sparse Fourier image (I5j . The 
crossover point between full and QTT versions of the modified Bini's algorithm depends 
on m and <x and can be even smaller, e.g. log 2 n ~ 1 7 for m small in modulus. 

5.2. Accuracy test for constant forcing 

We consider a simple problem for which the analytical solution is available, namely 
the one with constant forcing term. 



DJy(t)=my(t)+A, y(0)=y . 
The analytical solution is written in the following form 

y(t) = u F a (mt a ) + — E K (mt a ) - — , 

m m 



(21) 



(22) 



where E K is the Mittag-Lefler function H28ll27|| , which can be expressed and computed by 
certain (sometimes slow-converging) series. 

The accuracy verification results are shown on Fig. [6j We see that as the problem size 
grows, the accuracy improves until certain point and then the error start growing. This 
is explained by the machine threshold errors amplified by the condition number of the 
matrix A from ^ which is unbounded as n grows to the infinity. 



5.3. Accuracy of the Laplace transform 
Consider the following test equation 

D?y(t)=my(t; 



3 



y(0) = i. 



(23) 



Since the forcing term f (t) = t 3,/4 does not have a short Taylor series representation, this 
problem could be difficult for methods based on it. Unlike the previous example, the an- 
alytical solution in space domain is not available. Instead we can solve the problem using 
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Figure 7: The Laplace transform of the solution p5| l and its discrete approximations. Fractional parameter 
a = 0.8, mass m = — 1 . Left: fixed number of grid points n = 2 22 , different maximum time T. Right: fixed 
maximum time T = 1 3 , different Ti- 
the Laplace transform^ The Laplace transform of a function f (t) with the appropriate 
speed of decay is defined by 



F(s) = £{f(t)} 



e" st f(t)dt. 



(24) 



The Laplace transform of a convolution is the product of Laplace transforms, 



(f*g)(t) 



f(x)g(t-T)dT & £{f*g}(s) =£{f}(s)£{g}(s). 



This allows to simplify the equation (23) and find the Laplace transform of the solution, 



Y(s) 



1 



+ 



r(i.75) 



s 1 a (s<x — ml s h75 (s a — m) 



(25) 



The inverse Laplace transform y(i) = £ _1 {Y( S )} is given by the complex contour in- 
tegral and is difficult for numerical computation. However, we can easily compute the 
Laplace transform of the discrete solution £ y = Y in points {s k } using a rectangle quadra- 
ture rule, 



Y(s k ) « Y k = hY_ e-^yftj), tj = jh. 



(26) 



Then we compare Y(s k ) and Y k to establish the accuracy of the discrete solution y (tj) = y j . 
It is easy to see that equation ( |26] > contains three sources of errors, i.e. the ones of the 
discrete solution, of the quadrature rule and of the truncation of indefinite integral (24) 
to the finite interval [0 : T], T = nh. To compute Y(s) accurately for small s we should 
take T > s 1 log e 1 , where e is a machine precision error. To keep the quadrature rule 
error small, we should also use grids with small time step H. This is shown on Fig. [7j 
where the exact Laplace transform (25) is compared with its discrete approximations for 
different T and n. These factors motivate the use of very large grid size n and hence the 



4 History of the Laplace transform and other essential details can be found in, eg. [3 J. 



17 



log, |?(s k ) - t t \ log, IFOiJ - Fi| 




Figure 8: Accuracy of the Laplace transform | [25) given by discrete approximation |26j. Fractional parameter 
a = 0.8, mass m = — 1 . Left: fixed number of grid points n = 2 28 , different maximum time T. Right: fixed 
maximum time T = 10 7 , different n. 



QTT approach. It should be noted that the discrete Laplace transform ( |26] > is computed 
perfectly in the QTT format since the QTT-ranks of the exponent are all ones. 

Finally, on Fig. [8] we show the accuracy of the Laplace transform of the solution (25 1 for 
1 3 < s < 1 0° and for different T and n. It is clear that large problem size is essential for 
the accurate representation of the solution in the Laplace transform space. 



6. Conclusions and future work 

We present the new family of algorithms for the solution of linear fractional ODEs. Our 
approach develops the framework of matrix algorithms for fractional calculus [41 J by em- 
bedding the QTT tensor decomposition inside matrices, as proposed in [33J. The proposed 
algorithms works on matrix level and can be formally applied for the inversion of any tri- 
angular Toeplitz matrix, as well as the one obtained by discretisation of a linear fractional 
calculus problem. The workhorse of the inversion algorithms is the discrete convolution 
and/or Fourier transform of vectors given /approximated in the compressed QTT form. 
The success of the proposed algorithms, however, is based on the representability of the 
initial matrix and intermediate vectors arising in computations in the QTT format with a 
modest accuracy. 

As the motivating example we consider a simple linear fractional differential equation 
which reduces to the weakly singular convolutional Volterra equation with the Abel-type 
kernel. The QTT approximation method benefits from both the smoothness and decay 
of the Abel kernel, which results in efficient QTT-representation of problem matrix with 
the accuracy up to the machine precision. As shown by numerical experiments, the QTT- 
ranks of the intermediate vectors in the proposed algorithms remain bounded or grow 
slowly with the problem size. As the result, our algorithms of the inversion of triangular 
Toeplitz matrices demonstrate sublinear o(u) complexity, which falls down to the com- 
plexity O(log 2 n) of the superfast Fourier transform in certain cases. For our implementa- 
tion the crossover point with the standard algorithms based on the FFTW library for the 
considered experiments is 1 7 < log 2 n < 21 , i.e., the developed methods give not only the 
asymptotical benefit, but also a practical speedup for the problems of moderate size. 

The proposed approach opens a new class of algorithms for the fractional calculus, 
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i.e., methods of sublinear complexity. The developed techniques can be applied to the 
fractional equations with several differential operators of different order. They also can be 
generalised to fractional PDEs in two and more dimensions and to the nonlinear fractional 
problem. This would be the topic of further work, which will be reported elsewhere. 
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A. Fractional differential operators 

We begin by presenting established definitions of the fractional Riemann-Louiville op- 
erator, the fractional Riemann-Louiville derivative and a modified form of the fractional 
derivative — the Caputo derivative. These definitions can be found in a variety of ealier 
sources, including [12 j and Q40]| . 

Definition 1. /I12I/ Let a e R + . The operator J* defined on Li [a, b] by 

(t-s) a - 1 g(s)ds 

for a < t < b, is called the Riemann-Louville fractional integral operator of order a. For a = 0, 
we set J°L = I, the identity operator. 

Definition 2. /I12I/ Let ot 6 IR + and p = [a] . T7ze operator D*, defined by 

D«g := DT Q -*g 

z's caZ/erf zTze Riemann-Louiville fractional differential operator of order oc. For oc = 0,we set D° := I, 
f/ie identity operator. 

Definition 3. /I12I/ Assume that a > and zTzai g z's swcfo tfzaf D* (g — T p _! [g; a]) exists, where 
p = [a] and T p _ 1 [g; a] z's the Taylor polynomial of degree p — 1 for the function g about the point 
t = a; T p _i [g; a] := Ofor p = 0. Then we define the function D" Q g by 

Df Q g :=D:(g-T H [g;a]). 

The operator D* a z's called the Caputo differential operator of order a. 

We have chosen a problem defined in terms of the Caputo derivative because it allows 
us to specify non-homogeneous initial conditions for our test equation and thus it is more 
advantageous for modelling real- world phenomena ||12||. This allows us to draw compar- 
isons with existing work (such as [T¥| and ||T7|) which focus on the Caputo form of the 
derivative. 

We also note that for our range of a we will always have p = [a] = 1 . Also, as is 
sometimes the case in the literature we will omit our starting value a = for t from our 
notation; i.e. we will write J" for J", D a for and for D" Q . 

Solutions to fractional problems such as |l]) are formulated as functions of the Mittag- 
Lefler function (first defined by Mittag-Leffler in ||28|, [27]). A discussion of the relevant 
properties may be found in | jT2| and H40] ]: we re-present the details necessary to this paper 
below. 

Definition 4. Let a > 0. The function E a defined by 

OO j 

E«(z) := Y_ w . — — pr 
whenever the series converges is called the Mittag-Leffler function of order a. 



J?9(t) := f 
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Definition 5. Let <xi,<x 2 > 0. The function E ai )Ct2 defined by 



Eai ,a2 (^) 



j=0 



r(jai + a 2 ) 



whenever the series converges is called the two-parameter Mittag-Lejfler function with parameters 
ot] and 0L2. 

Note that E K (z) = E«,i(z). 

Theorem A.l. (/O) Consider the two-parameter Mittag-Leffler function E ai , a2 /or some cxi , ct 2 > 
0. TTie power series defining E ai )K2 (z) zs convergent for all z G C. 

Theorem A.2. ( 71121/ ) Lef a > and AeR. Moreover define 

y(t) :=E a (At a ), x > 0. 

TTzen 

D«y(t)=Ay(t). 

It is straightforward to apply the more general theory and results in texts such as | |12 || 
(which have their origins in analogous results from classical calculus) in order to prove 
that a solution exists for over a finite range of t. We present such adapted results. 

Theorem A.3. Let K > 0, b > 0. Define G := {(t,y) : t G [0, b], |y — yd < K\and let z : G — > 
M, swcfo £/zaf z(t,y(t)) = my (t) + f(t) under aboue conditions, be continuous. Furthermore 
define M := sup (ty)gG |z(t,y(t))| and 



B = 



min < b, 



KT(a+l) \ « 
M 



M = 
otherwise 



TTien fere ex/sis a unique function y G C[0, B] solving Q. 



In order to prove Theorem A.3 we first need to prove the following lemma: 



Lemma A.4. Assuming the conditions of Theorem A.3. the function y G C[0, B] is a solution of 
the initial value problem ([!]) if and only if it is a solution of the Volterra integral equation 



y(t) =yo + 



Via.) 



iot-1 



zfslds. 



(A.1) 



Proof of Lemma A.4: Assume y(t) is a solution of (A.l I. Writing (A.l I in operator form 
we have 

y(t)=y + J a z(t) (A.2) 
and subsequently applying D* to both sides yields 

D?y(t) = D,"yo + D?J"z(t) 
= + z(t,y(t)). 
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Thus, recalling that z(t,y (t)) = my(t) + f (t), y(t) must also solve Q. Proving the condi- 
tion in the other direction, we now assume that y(t) is a solution to ( |A.l >. Recalling that 
z{i] G C[0, B] we may write ^ as 



z(t,y(t)) 

Applying J to both sides yields 

Jz(t,y(t)) 

Now applying D 1 a to both sides we have 

D 1 - a Jz(t,y(t))=y(t)-u , 
which may be rearranged to give the precise Volterra equation we require: 

y + J a z(t,y(t))=y(t).D 



D«y(t) 

D a (y-y ) (t) 
DJ 1 - a y(t)-DJ 1 - a y 



JDJ 1 - a y(t)-JDJ 1 - a y 

r^yw-r^yo 



Proof of Theorem [A3) Suppose that M = 0; then z(t,y(t)) = 0V(t,y) G G. For this 
case, y : [0, B] — > M such that y (t) = y is a solution of the initial value problem ([T]). 

Otherwise, supposing M ^ 0, we must use our Lemma A.4 which asserts that ([T] is 
equivalent to the Volterra equation (A.l I. We introduce the set U := {y G C[0, B] : ||y — yolloo 
U is a closed, convex subset of the Banach space of all continuous functions on [0, B], 
equipped with the Chebyshev norm. Hence, U is a Banach space also. Since y = y G U 
also, we conclude that U is non-empty. We define an operator R on U by 



<K}. 



(Ry)(t) :=y + 



1 



r(«) 



t-sr-^yts^ds. 



We can now rewrite our Volterra equation ( A.l[ > as 

y = Ry 



(A.3) 



(A.4) 



and our task of proving the existence of a solution to (A.l I (which is equivalent to proving 
the existence of a solution to our original problem ([l I) now becomes one of proving that 
the operator R has a fixed point. Consider that, for < ti < tz < B, we have 



|(Ry)(ti)-(Ry)(t 2 )| 



1 



< 



< 



r(«) 

M 



r*2 



\a-1 I 



(t 2 -s) a - 1 z(s,y(s))ds 
r'2 



(since a < 1 
2M 



as tz — > ti . 



a 

(t 2 -tO a 



- (t 2 - s 
1 < => (ti - s 



ds + 

a-l 



(t 2 -s) a - 1 ds 



ti 

> (t 2 
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In addition, for y G U, t G [0, B] we have 



|(Ry)(t)-y | = 



< 



Via) 



MKr(a+ 1 



r(a+ i)M 



= K. 



Thus, y G U => Ru G U. In order to show that we have a fixed point we must now show 
that R(U) := {R(u) : u G Y} is a relatively compact set. For w G R(U) we have, for all 
TG[0,B], 



|w(t)| 



l(Ry)(t)l 
1 



If |t 2 - ti | < 6 then 



r(a) 
< -yo + K. 



|(Ry)(ti)-(Ry)(t 2 )|< 



ia-1 



z(s,y(s))|ds 



2M6 a 
r(a+V 



thus the set R(U) is equicontinuous (due to the right hand side's indpendence of u, t 2 , t] ). 
We can therefore apply the Arzela-Ascoli Theorem [U to conclude that A (U) is a relatively 
compact set and consequently apply Schauder's Fixed Point Theorem | |44 || to conclude that 
R has a fixed point. By construction, our fixed point solves the original problem ([I]). All 
that remains is to prove that our solution (whose existence we have just proved) is unique. 

Suppose we have a second solution (i.e. a second fixed point), y. If we consider 
\\V ~ V I loo t nen repeated use of Ru = y and Ry = y yields 



< 



(|m|B 



ro + <xj; 



\y-y\ 



To confirm uniqueness then, we just need to confirm convergence of the above; and 
what we have on the right hand side of the inequality is the power series definition of the 
Mittag-Leffler function E K (|m|B a ), which we know converges, by theorem A.l Our proof 
is therefore concluded. □ 
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